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ABSTRACT 

bo" 

O ! We study the evolution of an accelerating hyperrelativistic shock under the 



O 



presence of upstream inhomogeneities wrinkling the discontinuity surface. The 
investigation is conducted by means of numerical simulations using the PLUTO 
code for astrophysical fluid dynamics. The reliability and robustness of the code 
are demonstrated against well known results coming from the linear perturba- 
tion theory. We then follow the nonlinear evolution of two classes of perturbing 
upstream atmospheres and conclude that no lasting wrinkle can be preserved 
indefinitely by the flow. Finally we derive analytically a description of the geo- 

^ ' metrical effects of a turbulent upstream ambient on the discontinuity surface. 
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Subject headings: hydrodynamics - shock waves - corrugation instability 



Introduction 



r^ ■ There seems to be strong evidences that Gamma-Ray Bursts (GRBs - see Piran (2005) 

^yQ • for a review) involve flows of dense shells thrown by a dying compact star in the ambient 

O ■ medium at Lorentz factors P > 10^ — 10^. When the ejecta impact on the surrounding 

matter carried by a pre-existing stellar wind, a shock is formed and begins to propagate into 
^ ' a decreasing atmosphere, circumstance which leads, for sufficiently steep density profiles, to 

c^ ■ the shock acceleration. The length scale fcg on which the stellar atmosphere rarefies maybe 

reasonably much smaller than the distance from the center of the star, thus justifying the 
approximation of planar symmetry in studying the shock evolution. 
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The problem of a 
both in its Newtonian 



jlast wave moving; into a decreasing atmosphere has b een analyzed 



Grover fc Hardv 



1966 



Gandel'man &: Frank-Kamenetskii 



Haves 



1968) and relativistic (Blandford &: McKee 



1956 



Sakurai 



Perna fc VietrJbooilNa^avama fc ShigevamalKsl i lPan fc Sarill2006l : ISari|[2006h 



960l:lRaizei]ll964J: 



1976; 



Best &: Sari 



2000 



regimes, 



and several self-similar solutions have been found for the flow in both power-law and expo- 
nentially shaped density profiles. Despite the importance of the issue, very few papers 
have been spent t o study th e stabi lity o f the system subject to wrinkling perturbations. In 
Newtonian regime IChevalierl (Il990l ) and ILuo fc Chevalieii (Il994r) studie d the exponential at - 
mosphere; power-law profiles have been considered in lSari et al\ (120001 ). IWang et al\ ( 120031 ). 
while taking into account relativistic effects, were not able to find any self-similarity in the 
perturbation analysis of a spherical blast wave propagating in a power-law atmosphere. 



Palma &: Vietril (120061 ) performed a linear stability analysis of a highly relativistic planar 



shock propagating in an exponential atmosphere and retrieved a self-similar solution for the 
first order problem. They obtained that, at least in the small perturbation limit, with respect 
to what happens in the Newtonian regime, the corrugation wavelength k~^ can drop by a 
factor of F still giving rise to no sensible restoring effect in the flow, a behaviour reminiscent 
of the infinite wavelength case, even for small ratios ko/k. This allows the instability of the 
downstream energy density to persist, thus delaying the saturation phase. 

Of course, as the instability exits the linear regime, a numerical approach has to be 
adopted since the arguments that exclude the arising of stabilizing phenomena in the flow 
may become weaker and not so pertinent. 

The plan of the paper is as follows. In ^ we rev iew t he analytical propertie s of the 
problem extensively discussed in iPerna fc Vietril (120021 ) and iPalma fc Vietril (120061 ) : these 
predictions provide useful benchmarks to which our numerical scheme can be compared. In 
§31 we describe the code used to integrate the relativistic hydrodynamics (RHD) equations; a 
specific subsection is reserved to explain how we overcome the relevant technical difficulties, 
enlists all the major tests through which we run the code before declaring it reliable for our 
purposes. In ^ we tackle the central subject of the paper, thus reporting several results of 
simulations dealing with nonlinear variations of the perturbations induced into the system 
in §H Lastly, in §H1 we will show that self-similarity is reached fast enough to allow for 
an analytical expression for the shock speed in a quite arbitrarily shaped atmosphere. By 
means of such a result we will develop a technique to calculate the shock position without 
making time expensive simulations and will apply it to describe the shape evolution of a 
planar shock impacting a turbulent upstream. Conclusions are drawn in the shape of an 
excursus in ^ 
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Self-similar solution 



In this section we summarize the predictions, analytically derived in Perna & Vietri 
(2002) and in Palma & Vietri (2006), whose accuracy we will check in the following. In the 
first one the self-similar solution for an accelerating hyperrelativistic shock propagating in a 
planar exponential atmosphere is derived. Assuming an ambient density given by 

p{x) = p,e-^'\ (1) 

dimensional and covariance arguments impose the self-similar shock speed V{t) (hereafter 
we will pose c = 1) to satisfy 

-a 2 °^(l-\/)(l + K)) V^Vo' ^' 

where Vq is the shock speed at time t = and a is a dimensionless (negative) constant to be 
determined by imposing a smooth passage of the flow through a critical point. As the shock 
enters the hyperrelativistic regime, Eqn. [2] becomes 

r(i)«r,expMi^«r.(^y'°. (3) 

-" \PiJ 

Here F is the shock Lorentz factor and the subscript i refers to the initial condition. 

In order to determine the value of a and the downstream profiles of the relevant hy- 
drodynamical quantities, the exact adiabatic fluid flow equations as well as Taub's jump 
conditions across the shock are considered in their highly relativistic limit. Chosen the 
self-similarity variable 

^ = ko[x~X{t)]T'{t) (4) 

(X(t) being the shock position), the hydrodynamics equations can be cast into self-similar 
form by means of the following separations of variables: 

72(x, t) = giOr\t) , e(x, t) = goi?(Or'+" W , (5) 

n{x,t) = zoN{0'r^+"{t), (6) 

with 7, e and n being, respectively, fluid local Lorentz factor, proper energy density and 
baryon number density (the first and last ones as seen from the upstream frame), go = Po/^f 
and Zq = tiq/V". Taub's jump conditions are satisfied simply by fixing 

^(0) = 1, i?(0) = 2, iV(0)=2. (7) 
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Solving the equations with respect to g{^), R{0 ^"^^ ^iO o'^^ finds that self-similar quan- 
tities satisfy the following Cauchy problem: 

^, ^ 2g [-2a(4 + a) + {2 + a){a - A^g] R 



a^ + {a- AO9 [-4a + (« - 4^^] 
g^4{a-40g-Ua-3a^] 



(9) 



" a^ + {a- 40g [-4a + {a - A^g] ' 

M2±a)Mz^. (10) 

Demanding the simultaneous vanishing of the numerators and denominators of Eqns. [8] and 
[9] at a critical point (thus specifying the "second type" nature of this self-similar problem) 
it is possible to find 

a = -(2 + 4/^3). (11) 



Lying on the zeroth order solution hitherto reviewed, iPalma Sz Vietril (120061 ) performed a 
linear stability analysis with respect to a shock wrinkle of wave number k. In the k/{kor) ^ 
1 limit, it is shown that causal phenomena transverse to the shock direction of motion 
cannot carry disturbances too far. This justifies the approximation of infinite wavelength and 
independent (zeroth order) evolution of each fiow column with slightly perturbed constants 
in the equation for the shock location. Strictly speaking, Eqn. [3] can be integrated to give 



If we perturb Ci we obtain 



(5Xocr°, (13) 

6e ex /2'(0r^+"(t) , (14) 

Sn oc iV'(Or^+"(t) , (15) 

57^ oc g'm\t) . (16) 



Alternatively, perturbing r,: 



6X = 6ci oc r-2 , (17) 

5e ex (4i?(0 + A^R'iO - c^R'iO) r'+"(^) , (18) 

6n oc (4iV(0 + 4eiV'(0 - aN'{0) r'+"(t) , (19) 

5^ (X (2^(0 + 2^g'{0 - f ^'(O) r^(t) . (20) 

It is clear that the first mode is the most severe and thus physically relevant for the instability. 
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Nevertheless they perform the full perturbation analysis which takes explicitly into 
account the transverse mixing between adjacent columns, thus also obtaining a complete 
description for the y-component of the four velocity: 



5uy (X gym'-\t) 



(21) 



with s being the parameter which selects the strong {s = 3) or the weak (s = 1) mode and 
Qy the self-similar profile satisfying 



g'gy [aR' + AgR{s - 3)] gy - ia{k/ko)^Ri 



2g 



9yW 



2R[{a-A0g-a] 
ik 



V2ko 



(22) 



(23) 



Here gy is purely imaginary since it is 7r/2-shifted with respect to the other perturbations. 

In the following two sections we will try to numerically recover near all the theoretical 
results stated above. 



3. Numerical Setup 



Numerical simulations are carried out by solving the equations of number density and 
momentum-energy conservation, i.e. 



dn -L 



[nv] 



0, 



Q'ffl -» 

— - + V ■ imv + p) 
at 



dE 
'dt 



+ V- m = 0, 



(24) 

(25) 
(26) 



where v is the fluid velocity, m = nmphVv and E = nmphV — p are, respectively, the 
momentum and energy density {rrip being the proton mass). 



Proper closure of Eqns. [24l^[26]is specified in the form of an equation of state (EoS), 
relating the specific enthalpy h = 1 + e + pT / {nrrip) with pressure p and (specific) internal 
ener gy of the flui d e. For a relativistic perfect fluid, the desired closure is given by the Synge 
gas (ISyngd 119571 ). For a single-specie fluid given by a mixture of protons and electrons, 
th e equation of state c an be approximated by an analytical expression late ly presented 



m 



Mignone et al\ (120051 ) and further discussed in iMignone fc McKinneyl (120071 ) : 



h 



9 + ^/^02 + 1 



(27) 



-6- 



where = p/{nmp) is a temperature-like variable. Compared to the ideal gas EoS with 
constant adiabatic index Tg for which the enthalpy takes the form h = 1 + Tg/lTg — 1)0, 
Eq. fl7r|) yields the correct asymptotic limits for very high (0 -^ oo) and low (0 —>■ 0) 
temperatures, r educin g to an ideal EoS with Tg = 4/3 and Tg = 5/3, respectively. In 



Mignone et al\ ( l2005l ) it is shown that this expression differs by less than 4% from the 
theoretical prescription given by the Synge gas. Since the equation of state is frequently 
invoked in the process of obtaining the numerical solution, computational efficiency issues 
largely entitle to the use of an approximated relation. 

The conservatio n laws (Eqns. [24l-[26l) are solved using the relativistic module available 



in the PLUTO code (JMignone et a/.ll2007l ). PLUTO is a Godunov-type code offering a variety 



of computational strategies for the numerical solution of hyperbolic conserva tion laws in one, 



two o r three dimensions. For an extensive review of such techniques see (JMarti &: Miiller 



20031 ) and references therein. Being Riemann-solver based, it is particularly fit for the 



simulation of high-mach number fiows, as it is the case here. For the present appl ication. 



we em ploy second-order accuracy in time by using characteristic backtracing (see IColella 



( 119901 )) and linear interpolation with 2"^^^ order limited slopes. This scheme yields a one-step 
time integration by providing time-centered fiuxes at zone boundaries, computed by solving 
a Riemann problem with suitable time-centered left and right states. For the one- and two- 
di mensional simulat i ons pr esented below, we adopt the approximate HLLC Riemann solver 



of iMignone fc Bodol tOOH ). 



3.1. The Choice of the Reference Frame 

Before presenting our numerical results, we discuss how we faced a number of numerical 
issues. 

Let us begin by considering the fate of an upstream slab one length scale long: due to 
the highly relativistic shock compression, it will be roughly resized by a factor F^. In order 
to justify the hyperrelativistic approximations assumed above, we would be willing to deal 
with F at least as big as 10; even higher Lorentz factor are involved with realistic models 
of GRBs, wherein the compactness problem solution imposes F to be one or two orders of 
magnitude higher. 

However, such large Lorentz factors demand an increasingly high resolution, if one 
wishes to properly capture dynamics of the slab profile one it enter. This requirement 
becomes even more severe if the evolution of the perturbations has to be followed accurately. 
In this case, in order to overcome spurious numerical fiuctuations, a resolution of thousands 
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computational zones per length scale is needed. From these considerations, we conclude 
that adopting a static uniform grid would result in extremely inefficient calculations. To 
overcome this limitation, a reasonable alternative is to resort to adaptive mesh refinement 
(AMR) techniques, thus providing adequate resolution on the regions of interest. Even in 
this case, however, we still have to face a subtler problem. 

It is known that relativistic shock- capturing codes may suffer from excessive dissipation 
when a region of fluid with exceedingly large inertia interact with a stationary fluid adjacent 



to it (JMignone et a/.ll2005l ). This is actually the unfavorable situation we are coping with 
since, in the upstream rest frame (URF from now on), an ultrarelativistic shock advances 
in a cold, pressure-less static gas. In this reference frame the jumps of the hydrodynamical 
variables (in particular the energy density) across the front are maximized, leading to an 
excessive smearing of the shock proflle. The deficiency is inherent to any finite difference 
method attempting to solve the fluid equations on meshes of flnite width. Indeed, even a 
flrst-order upwind discretization of the scalar advection equation with linear constant velocity 
c> 0, 

_^_ 3_ ^ ^ 1 1 1 ^ Q 28 

At Ax ^ ^ 

shows that u satisfles exactly another convection-diffusion problem, namely 

where the term in brackets on the right hand side must be positive for stability issues. Thus, 
roughly speaking, the magnitude of the second derivative provides a rough estimate of the 
diffusion introduced by the numerical algorithm. The situation does not improve with the 
employment of higher order methods, since the accuracy reverts to flrst order in proximity 
of a discontinuity anyway. This conclusion is supported by several numerical experiments 
(not shown here) showing that the largest dissipation terms, taken to be proportional to the 
magnitude of the second derivative of the hydrodynamic variables, result in the frame of the 
upstream fluid whereas are minimized in the shock frame. 

In this respect, it is thinkable to study the shock evolution in its initial instantaneously 
comoving frame (IICF). In fact, in such an inertial reference frame, shock compression results 
to have a modest factor of 3 (i.e. an upstream slab of unitary length will be resized by a 
factor 3). Moreover, in the same frame, due to the favorable Lorentz factor composition law, 

r' = rro(i + /5/5o), (30) 

the shock will hardly become hyperrelativistic even in the late acceleration stages. This 
allows to follow the long-term evolution of the downstream self-similar lengths as well as the 
upstream length scales with a comparable number of points. 



The price one pays for reducing in such a drastic way the computational cost consists in 
a quite complex procedure to recover a snapshot of the system as seen by an observer at rest 
with respect to the upstream. Due to the relativistic non-absoluteness of simultaneity, we 
had to make a collage with several (ideally infinitely many) pieces of IICF snapshots, each 
depicting a particular (ideally infinitely narrow) slab (normal to the x-axis) of the flow at a 
particular IICF time. 

In particular, referring to IICF quantities by means of primes, we chose as initial con- 
dition a planar shock located at Xq = Xq = moving rightward in a grid covering a unitary 
IICF length along the x-axis. We assumed a uniform downstream flow connected to im- 
mediately pre-shock upstream by usual Taub's jump conditions: such a choice corresponds 
to a shock which at t, t' < propagates in a uniform, cold (hence stationary) atmosphere 
which, at X = x' = 0, turns exponential. If the simulation lasts t'^ and we are interested in 
x' > x'l region (corresponding, in the URF, to a semi-infinite patch which closely follows on 
the left the hyperrelativistic motion of the shock) we can only inquire into shock evolution 
up to t < ro(tg + Vqx^). Let us consider now an array X' = {x'j} containing the distinct 
x-component of the grid. Lorentz transformed array of X' in the URF, 

X = ^+(3ot, (31) 

J- 

contains the x-components of the leftmost sector of the patch introduced above. In order to 
derive hydrodynamical quantities Q as measured in the URF at point Xi and time t one has 
to analyze the snapshot taken in the IICF at time t' = ro(t — PoXi)'- 

Q{x„t) = FiQ'{x'„t')). (32) 

Here F denotes the map functions transforming density, velocity and pressure from the IICF 
to URF. Since the time-marching algorithm evolves by discrete time steps, we performed 
a linear interpolation between the two set of quantities obtained by replacing t' in Eqn. [32] 
with, respectively, f""^ and t'" such that t'""^ < t'{t,Xi) < t'". 

Such a discussion can be easily extended to an AMR structure, having care to perform 
temporal interpolation between the finest level step times available at each spatial position. 



3.2. Simulation settings 

In all the simulations described in §1] and ^ the initial shock Lorentz factor Fq is set 
to 50. Similar results have been obtained by studying the evolution of shocks with different 
highly relativistic initial Lorentz factors. The upstream density of the atmosphere swept 
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up by the shock spans over two orders of magnitude, with a value of k^'^ ^11. The only 
exception concerning To is the simulation described in §4.21 and illustrated in Fig. [2l while 
in §5.21 we followed the shock evolution over ^ 7 length scales. 

In all the ID simulations described in this paper, an effective resolution of 2.56 x 10^ 
grid points has been reached by means of 8 refinement levels on a base grid with 10^ cells; 
the domain box was [0, 1]. 

The 2D simulations were performed on static grids with spatial resolutions indicated in 
Tabled 
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Table 1. Physical parameters and resolution adopted in the two-dimensional 

simulations. 



k/ko 


e 


rko 


[xb,Xe] X [yb,ye] 


Resolution 


Fig. 


4.8 


0.5 


- 


[0,1] X [-5,5] 


2 - 10^ X 5 - 10^ 


m 


24 


2 


- 


[0,l]x[-l,l] 


1.2 -10^ X 2.4- 103 


Eiiai 


4.8-102 


0.5 


- 


[0,1] X [-5-10-2,5-10-2] 


4 - 103 X 4 - 102 


Eiiioiiiiiin] 


4.8-102 


2 


- 


[0,1] X [-5-10-2,5-10-2] 


4 - 103 X 4 - 102 


[I1II1II3II1] 


4.8 -10^ 


0.5 


- 


[0,1] X [-5-10-3,5-10-3] 


10* X 4 - 102 


mm 


- 


3 


2.4- 10-1 


[0,1.6] X [0,2.285714] 


1.92- 103 X 1.38-103 


mmm 



Note. — The domain box is defined by the lower and upper coordinates [a;fc,a::e] (in the x 
direction), [yb,ye\ (in the y direction). Both the domain box and the resolution refer to the base 
computational grid used in the code. 
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Code Verification 



In the following subsections we will try to recover near all the theoretical results dis- 
cussed in §2. Such an "exercise" will provide us a powerful tool to test the code and the 
analysis procedure itself together with a measure of the reliability of those simulations having 
no clear theoretical counterpart (nonlinear perturbation, transient to self-similarity). 



4.1. Zero-th order solution: hyperrelativistic regime 

Let us consider the simulation of a shock propagating till a fixed point in a homogeneous 
atmosphere. As it enters a region with an exponentially decreasing density profile, the shock 
exhibits some inertia and its speed setups to the self similar value (Eq. Eqn. [3]) after a little 
while, see §6. 

Incidentally, we would like to point out that such a problem is totally scalable with 
respect to the length scale kQ^: having set c = 1, we are still free to choose the space (and, 
therefore, time) measurement unit. As a consequence, all the results we will obtain in this 
paper are completely independent of the specific value assumed by fco- 

As the shock advances into the stratified atmosphere, it appears possible to study both 
the spatial profile of the downstream hydrodynamical quantities and their temporal evolu- 
tion. 

The snapshot in Fig. [1] shows the density, pressure and Lorentz factor (normalized to 
immediate post shock values of 2, 2 and 1/2) in a small region immediately behind the shock 
front. The size of this region is a factor of ~ 2 smaller than the length traversed by the 
point originally marking the atmosphere change from homogeneous to exponential. We also 
overplot the curves obtained by direct integration of Eqns. |H]-[ini As clear from the plot, 
we obtain an excellent agreement. 
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Fig. 1. — The figure shows, from top to bottom, the spatial dependence of the dimensionless density n(^), 
pressure e(^) and squared Lorentz factor .g(^): theory predictions (soHd hnes) are compared with numerical 
results (40 crosses, stars and diamonds sample in the figure the numerical data). 
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On the other hand, one can check theoretical rules about Lorentz factor growth as 
function of time (or, equivalently, under the number of length scale swept up by the shock 
throughout the simulation). This can be done by means of several, sometimes equivalent, 
ways. Here we report only two of the most direct methods to be implemented (those our 
experience suggests to be likely the most robust ones). 

A first consistency check can be made by comparing the theoretical value of F predicted 
by Eqn. Owith the one computed from our numerical simulations. The latter can be recovered 
by solving Taub's jump condition with respect to the shock Lorentz factor once pre- and 
post-shock values have been identified. 

As an alternative, one can consider the value C,l of the self-similar variable corresponding, 
in previous fits, to the leftmost point plotted in FigdJ Strictly speaking, ^^ must be calculated 
by numerical inversion, for example, of -R(C) at the point Rl measured as the leftmost 
theoretical prediction in plot [1] Once C,l is known, one can recover the shock Lorentz factor 
connected to the simulation by inverting Eqn. HI 



^"^ (33) 



XL-X 

A comparison with the usual theoretical value completes the test. 

Perhaps it is a point worth of remark the fact that the latter sounds a bit more stringent 
test than the former, since explicitly assigns a fundamental role to the spatial profile of the 
downstream in determining the shock speed evolution, thus allowing a more complete point 
of view on the issue. 

Both tests provide an excellent agreement with the related predictions: deviations from 
theoretical rules appears to be nothing but numerical noise and at almost any time t > 
remain below a small fraction, typically less than 1%. 

As a measure of the code reliabihty we report what simulations predict about the pa- 
rameter a. From Eqn. [3] 

log^ 
a = -^ ; (34) 

substituting simulation values and averaging on several snapshot times we obtain 

a = -4.3102... , (35) 

a value which differs from the correct one for less than 0.02%. 
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4.2. Zero-th order solution: transrelativistic regime 

In this subsection we focus on tlie only prediction we have as long as the shock is neither 
Newtonian nor hyperrelativistic: the rule about shock speed given by Eqn. [2l We consider a 
shock with an initial Lorentz factor To = 1.1 and follow its evolution for the crossing of ~ 4 
length scales. In Fig. [2] we plot F as a function of traversed length (in units of k^^), together 
with the exact self-similar prediction (Eqn. [2]) and the run as expected if the shock would 
have been highly relativistic. 

The excellent agreement we observe in this plot completes the thorough picture about 
the zero-th order problem. 
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Fig. 2. — Evolution of the shock Lorentz factor as obtained by simulation (stars) together with exact 
self-similar solution (solid line) and its hyperrelativistic approximation (dashed line). 
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4.3. First order solution 



In principle, the linear perturbation a nalysis of a planar shock is a fully 2D prob- 
lem ( IChevalierl Il99d : iPalma &: Vietril |2006| ). However, in order to approach the problem 
from a numerical point of view, it is convenient to take advantage of the infinite wavelength 
approxi mation, whose applica bility in the hyperrelativistic regime has been discussed in de- 
tail in (jPalma fc Vietril |2006| ). Such a scheme enables an almost full investigation of the 
physically relevant phenomena, still allowing a reduced numerical cost. 

The idea can be summarized as follows. Firstly we perform the usual ID simulation 
of a planar, unperturbed shock wave, as described in previous sections. Then we carry out 
a second ID run by perturbing the upstream region with an overdense (by a factor e) bar, 
limited in extension to a fraction A of length scale (hereafter A ^ 1.32); since also this last 
simulation is ID, the reader should imagine such a bar indefinitely extended perpendicularly 
to the shock speed. Moreover - in order to avoid spurious structures, here as well as in the 
following section - we joined smoothly the perturbing bar to the background up to the 4th 
derivative by means of the factor cos'^lirkolx — x) / A], x being the center of the bar (hereafter 



X 



O.TO/cq ), thus obtaining the following upstream density profile: 



U(x) 



Poe 



-kox 



1 + e cos 



nko{x — x) 

A 



^'^+i 



X 



H 



X 



A 

Wo 



X 



(36) 



where H is the Heaviside step function. Perturbations to hydro dynamical quantities are 
simply obtained by subtracting term-to-term the values of the second set from the first ones. 
No 2D simulations were needed and the original resolution along the x-axis of the zero-th 
order analysis has been maintained. 

Nonetheless, we warn the reader that any spurious oscillation that might marginally 
affect the zero-th order profiles, will result here in a severe noise, when trying to recover 
a variable value as difference between two slightly different quantities. This justifies such 
a high resolution which prevents us from performing directly a 2D simulation: only this 
way we can keep the noise down to a reasonable threshold. In order to fix this problem we 
smoothed the data by performing a local regression using weighted linear least squares with 
a 2nd degree polynomial. In this way we were allowed to recover the physically relevant 
development of the perturbation by smoothing away less than 5% noise. 

As in the first subsection, the work can be split up into two stages: the study of the 
spatial perturbation profile and the test of correct temporal growth. 

Results for the spatial dependence of perturbations of the hydrodynamical quantities 
(at a fixed time) are reported in Fig. [3l according to theory predictions and to our numerical 
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scheme. The agreement is acceptable, specially considering the way such numerical results 
are achieved. 
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Fig. 3. — The figure shows, from top to bottom, the spatial dependence of the perturbations to pressure 
(5e(^), density (5ri(^) and squared Lorentz factor (57^(^) normahzed to the immediate downstream value: 
theory predictions (solid lines) are compared with numerical results (40 stars, crosses and diamonds sample 
in the figure the numerical data). 
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Concerning the instability rate, it is possible to compare the temporal dependence of 
theoretical and numerical perturbations by studying how the simulated amplitude - normal- 
ized to the expected growth - diverges from unity. In Fig. H] the results of such a study are 
reported: only a minor gap of few percentage points appears, thus stressing once again the 
good suitability of our scheme even for the subtle task of studying wrinkle perturbations to 
an hyperrelativistic shock. 
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Fig. 4. — Temporal evolution of perturbation amplitude normalized to the expected growth. Q is defined 
as the ratio between simulated and expected perturbation growth; consequently, any gap between theory 
and numerical test should have as a counterpart a departure of Q from 1, commensurate with the gap. 
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4.4. Finite wavelength wrinkles 

Here we just report the results of a 2D simulation (obviously much less resolved in 
x-direction than the ID previous ones) dealing with an overdense upstream bar not uni- 
formly extended along the y-axes as in the previous subsection. Instead, we impose a finite 
wavelength (A = 2'7ik~^) sinusoidal profile (along with periodic y-boundary conditions), thus 
allowing in the following section a direct comparison with the effects of an analogous bar 
inducing nonlinear perturbation to the system. 

li k ^ 4.8/^0, the effect of a bar of amplitude e ~ 0.5 and extended A are reported in 
Fig. [5j we emphasize that the finiteness of the wavelength imply a non-zero y component of 
the velocity. 
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Fig. 5. — Hydrodynamical quantities in the linear perturbation regime {k « 4.8A:o; £ ~ 0.5; A « 1). 
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In particular it appears noteworthy that the downstream profiles of Uy in Fig. [5l al- 
though not sufficiently refined to test the code strictly speaking, are nevertheless completely 
consistent with theoretical predictions given by Eqns. [22] -[23] for the strong mode s = 3. 



5. Nonlinear perturbations 

Once the robustness and consistency of the scheme has been demonstrated, we are 
allowed to study 2D problems heavily involving completely new phenomena or, at least, 
processes neglected in the small perturbation regime. 

In the following we present firstly an idealized problem, aimed at inquiring whether or 
not the instability reaches any sensible saturation point: the weak sinusoidal bar discussed 
in the previous section will be replaced by a more substantial one. To follow we will show 
what happens if a dense cylindrical cloud hampers the downhill path of the shock. 



5.1. Sinusoidal bars 

Aiming to study the nonlinear phase of shock perturbations, we will impose here the 
sinusoidal profile of the perturbing upstream on the density logarithm {i.e. Il{x,y) = p{x) ■ 
-^QEsmkyj ra^j^g]; than on the density itself - as done previously. In this way we are allowed 
to use larger amplitude perturbations which often imply a contrast of several orders of 
magnitude between overdense regions and adjacent vacua. As usual, in all the following 
simulations, the sinusoidal bars have a width Ak^^. 
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Fig. 6. — The density logarithm (left) and the parallel 4-velocity Ux (right) for k ss 24/co and e « 2. Please 
note that each panel refers to a different time, according to the different phases discussed in the text. 
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Let us begin by considering a long wavelength density bar {k ~ 24ko) with an amplitude 
of about two orders of magnitude (e ~ 2). Such an upstream inhomogeneity induces a strong 
corrugation (see left panel in Fig. El depicting the barionic density in a late phase of the 
bar shocking process) and highly nonlinear perturbations to downstream flow. In Fig. El for 
instance, the right panel shows Ux just after the shock emerges from the bar. The high-speed 
blob we observe just behind the crest is the relic of the flow corresponding to the acceleration 
phase in the bar vacuum. Once the shock emerges from the low-density region, the impact 
on the unperturbed atmosphere produces the reverse shock which separates the fast blob 
from the main discontinuity. 

Moreover, similarly to the linear case discussed in §4.4, the shock tends to fill the 
valleys present in its profile simply by means of something like a potential flow of matter 
along the discontinuity normal: Fig. [8] shows that the matter flows from the crest to the 
valley. However, at least on the explored time scales, such a long perturbation wavelength 
prevents almost perfectly the gap between crest and valley from a quick damping that the 
mechanism above described would cause to higher wave-number wrinkles (see the sequence 
of snapshots in Fig. [7j). 
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Fig. 7. — The density logarithm for k « 24fco and £ « 2. 
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Fig. 8. — The transverse 4-velocity Uy for k « 24fco and e « 2. 
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At shorter wavelengths, the system evolves in a quite different way. 



-29- 



Log(|0), t = 5.000 



■ 


4.3 


■ 


3.6 


- - 


2.S 


- - 


2.1 


- - 


1.3 


■ 


0.6 


■ 


-0.2 





H 


4.3 


■ 


3.5 


- - 


2.8 


- - 


2.0 


- - 


1.3 


■ 


0.5 


■ 


-0 3 




Fig. 9. — The density logarithm for fc « 4.8 • lO^fco and e « 0.5. 
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Fig. 10. — The density logarithm for fc w 4.8 • 10^/co and e « 0.5. 
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Fig. 11. — The transverse 4-velocity Uy for fc w 4.8 • lO^fco and e « 0.5. 
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Fig. 12. — The transverse 4-velocity Uy for fc w 4.8 • lO^fco and e « 0.5. 
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Let us focus our attention on a bar of A; ^ 4.8 ■ lO^fco and e ~ 0.5. One can argue 
that, due to the reduced transverse distance between equally out of phase flow columns (or, 
equivalently, thanks to the higher gradients involved), the evolution observed in the previous 
run resembles the present case in slow-motion: the more fc/^o grows, the faster the evolution 
gets. In fact, looking at the sequences of snapshot in Figs. [9] - [TO] and [TT] - [T2] which depict 
the time evolution respectively of barionic density and My, it is possible to see that a situation 
similar to the last snapshots in both Figs. [7] and [8] here is reached on a reduced time scale. 
What we observe in all its progression is the sharpening of the valley, which starts being 
U-shaped and evolves in the shape of a V. At that point in practice we have two distinct 
shocks which go to intersect, with the resulting formation of two secondary shocks in the 
downstream. Such a X-shaped structure evolves with the secondary shocks which advance 
toward the adjacent crests. In this context, the valley closes the gap on the crest in the 
short and within another few fractions of length scale comes to overtake the leading front of 
the shock. At this point the play of the flow columns repeats with reversed roles. However 
it should be plain that, with each role reverse, two main things happen: first, a new layer 
is added to the existing pattern of hydrodynamical fluctuations the flow advects far in the 
downstream; second, the wrinkle amplitude gets smaller and smaller, thus coming to restore 
the original zero-th order solution. Having these facts in mind, it is possible to give an 
estimate of the time Tsm needed b y the shock, once it com es out of the perturbing region, to 



restore the original planar shape. iPalma &: Vietril (120061 ) most clearly discussed the scaling 



of such a time, so that we can say T^^ oc T /k. The coefficient of proportionality can be 
estimated from the simulation: if we say that planarity is restored in the last snapshot in 
Fig. [in] and remember that the shock comes out of the perturbing bar (of width Afeg) at 
t ~ 11, we obtain 

Ts^ ~ 25^ . (37) 

k 
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Fig. 13. — The density logarithm for fc « 4.8 • lO^/co and e « 2. 
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Fig. 14. — The density logarithm for fc w 4.8 • 10^/co and e w 2. 



-36- 



u» t - s«» 




0-04 
0.02 




■ 


124.2 
103.5 
82.8 


>. 0.00 


^^ ^^^^^^^^^^^^^^^^1 


- - 


62.1 


-0-02 


^^^^^^^^^^^^^^^^^^^^^^^1 




41.4 


-0.04 


^^^^^^^^^^1 


m 


20.7 
0.0 



Fig. 15. — The parallel 4-velocity Ux for fc w 4.8 • 10^/co and e w 2. 
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Fig. 16. — The parallel 4-velocity Ux for fc w 4.8 • 10^/co and e w 2. 
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We also include some snapshots of two other simulations showing close analogies with 
the previous one. The first and most spectacular one deals with k ~ 4.8 ■ lO^fco and e ~ 2, 
and presents clear evidences of Kelvin-Helmotz instabilities; compare Figs. [T3] - [IH and [TS] 
- Uni at early times, behind the valley, due to the stopping presence of the overdensity, a 
thin, slow layer courses the unperturbed, fast flow, giving rise to the instability. The second 
one, realized with k ^ 4.8 • lO^ko and e ~ 0.5, quickly comes to a restoring of the shock 
surface planarity (Fig. [T7]). The complexity of the fluctuation pattern described above can 
be seen from this last simulation: Fig. [TS] shows the tightly arranged warp of Uy already at 
an intermediate evolutionary phase. 
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Fig. 17. — The density logarithm for fc « 4.8 • lO^ko and e « 0.5. 
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Fig. 18. — The transverse 4-velocity Uy for k w 4.8 • 10"^fco and e ~ 0.5. 
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5.2. Cylindrical clouds 

The perturbations we are going to deal with here complement the ones we explored 
above. Indeed, if the sinusoidal bars provide a quite thorough description of those phenomena 
occurring during the sweeping up of a smoothly inhomogeneous upstream, a cylindrical 
overdensity can well represent the sharp contrast of a typical cloud in a clumpy circumburst 
medium. 

We present here the results of a simulation with a homogeneous cylindrical cloud (with 
axis parallel to the 2;-axis) of radius r ^ 2.4 ■ lO'^kQ^ placed in the upstream of the usual 
planar shock. The cloud density is larger by a factor 10^ than the unperturbed upstream 
one. 

In order to save computational time, we considered only half cloud and substituted 
periodic y-boundary conditions with reflective (rigid walls) ones. 

Figs, dni and [20] show the main evolution phases of the system: closely related to the 
first widest sinusoidal case, the unperturbed shock tends to fill the valley from the boundary 
of the cloud. However, in this case, because of the steeper wrinkle in the shock, the non- 
adiabaticity of the flow across the discontinuity surface allows a vortex to develop and to be 
advected downstream (Fig. [2T1) . The study of the vortex dynamics and its relevance with 
regard to GRBs physics will be discussed in a forthcoming paper. 
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Fig. 19. — The density logarithm for r « 2.4 • 10~'^kQ and e « 3 
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Fig. 20. — The parallel 4-velocity u^, for r « 2.4 • lO^^fc^^^ and e « 3 
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Fig. 21. — The vorticity V x {h-/v) for r w 2.4 • lO'^k^^ and e w 3. 
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6. Turbulent ambient density 

We discussed before the evolution of a shock in a self-similar regime under the effect of 
some perturbing agent. We now wish to investigate how fast such a self-similar regime is 
reached. 

Let us consider the usual shock with initial Lorentz factor Fq which, at t = 0, comes 
out of a homogeneous upstream and starts to propagate into an exponential atmosphere. In 
Figs. [22] -[21] we compare (for three different values of Fq spanning a wide range of relativistic 
regimes) the simulated Lorentz factor as a function of the length scale fraction covered by 
the shock with the prediction given by the self-similar theory (Eqn. [3]). 
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Fig. 22. — Shock Lorentz factor evolution as function of the distance traveled into the exponential atmo- 
sphere for To = 2 (stars for simulation, solid line for self-similar prediction, dashed line for the hyperrela- 
tivistic limit of the latter). 
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Fig. 23. — Shock Lorcntz factor evolution as function of the distance traveled into the exponential atmo- 
sphere for Fq = 10 (stars for simulation and solid line for self-similar prediction). 
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Fig. 24. — Shock Lorcntz factor evolution as function of the distance traveled into the exponential atmo- 
sphere for Fq = 50 (stars for simulation and solid line for self-similar prediction). 
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It is quite a significant fact that self-similarity is reached almost immediately after the 
exponential length scale switch (in our scheme, the length scale of the atmosphere switches 
from cxD for x < to kQ^ for a; > 0), with a very weak dependence on Fq: one senses that 
larger Lorentz factors help the shock to follow closer the self-similar run of the acceleration 
in the first 15 - 20% of length scale after the atmosphere transition. 

We now show that such a property allows us to consider the shock behavior as uniquely 
ruled by atmosphere density value at the point of interest. Obviously the II type nature of 
this self-similar problem will play a fundamental role here, since we are actually requiring 
the shock to have an infinite piston behind it which compensates for the indefinite energy 
supply in the upstream. 

Let us consider a generic upstream atmosphere density profile n(x). We assume n(a;) to 
be a monotonic decreasing function. It is possible to approximate n(x) with a finely broken 
line made up by a set of segments of exponentials. At least as self-similarity is reached 
on temporal - and thus spatial - scales smaller than the spatial scales required by n(x) 
to appreciably depart from the local tangent exponential, one is allowed to determine the 
infinitesimal increases of F by means of Eqn. (3) 

r.^.r(^y\ ,38) 

Integrating step by step Eqn. [38l the generalization to an arbitrary density profile of Eqn. [3] 
is easily obtained: 



r.r.^-j . (39) 

As a result, shock speed will depend only on the initial Lorentz factor and on the local 
density value. 

Having this fact on our mind, we can calculate how much a shock propagating in an 
atmosphere n(a;) will fall behind an initially identical one (also IIq = po) which propagates 
- unperturbed - in the usual exponential profile p{x). Calling Xi the former's position and 
Xo the latter's one and imposing Xo{t = 0) = Xi(t = 0) = 0: 

2kot 



V fr ^oXo\ ^ exp^ 

Xo = .^roexp-^^l-^^; (40) 



fJ^^ 



2/a 



X. ^ . I To (S^i) V" U I - ii|^ . (41) 



2rj 
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Subtracting Eqn. HO] from Eqn. UT], one obtains the equation for the delay: 



dt 



-AX 



2TI 



2knt 



exp 



a 



_Po_ 
U{t) 



2/a 



(42) 



Eqn. |42]is an interesting resuh, since it may prove to be useful for accurately estimat- 
ing the shock position in an arbitrary atmosphere without wasting any time in full blown 
simulations. 



As an example, let us write the density profile as 

n(x) = p{x) {1 + eSn{x)) . 



(43) 



If we are dealing with a slightly perturbed upstream, e will be much smaller than 1, thus 
justifying the following simplifications to Eqn. |42l 

e6U{t) exp ^ 



d 



AX '.. 

dt 2To 



1 exp^<|l-[l+£OT(t)]- 



a 



aTl 



(44) 



Let us suppose now that a planar shock encounters a turbulent upstream. According to the 
approximation of independent evolution of each flow cylinder we extensively discussed above, 
we can derive the statistical properties of the shock wrinkles at each x (or, equivalently, at 
each t). As a starting point we can imagine that each upstream cylinder perturbation 
(identified by a pair {y,z)) is a particular realization of a power spectrum A{k), such that 



A{k)= 6U{k,y,z) 



. Vy, z: 

5U{x,y,z) 



dkA{k) cos {kx + S{k)) 



(45) 



where the phase 6{k) is a random variable determining each realization with a probability 
distribution given by 

P(sik))^^imui^2i^m. (46) 



27r 



Integrating Eqn. 

AX{t) = 

and substituting Eqn. HS] we obtain 

"^ eA{k) 



e6U{T) exp ^dr 



AX{t) 



arg JO 



aTl 



cos (kr + S(k)) exp drdk . 

a 



(47) 



(48) 



Here two cuts have been introduced in order to exclude from the computation turbulence 
wave-length larger than a; ~ t itself (infrared cut k = t^^) or smaller than the scales reached 
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by transverse diffusive phenomena smoothing out high wave-number wrinkles (uhraviolet cut 
K = 25rt~^ - see Eqn. 1371) . From Eqn. HHJit is obvious that AX is itself a random variable 
given by the sum of (infinite) random variables, each of them identified by the parameter k 
and depending on the random phase 5{k). The central limit theorem completely characterizes 
(from a statistical point of view) AX as a random variable Gaussian distributed with an 
average over the ensemble of cylinders (AX(t)) = and a variance a = (AX^) given by the 
(infinite) sum of the variances da{k): 



rfa_r[^/o-s(fcr + 5)exp^rfr 



dk 2n 

whence, summing over k, it results: 



dS 

- , (49) 



K 



A^{k) [1 + exp ^ - 2 cos(A;t) exp ^] 



The characterization of the shock position distribution provided by Eqn. [50] may rep- 
resent a benchmark for several applications. Minor modifications should be sufficient, for 
instance, to obtain an estimate of the features in the fluctuations in the afterglow light 
curve, provided a model of upstream turbulence is given. Such a kind of study is considered 
of great interest, particularly in light of the recent Swift observations: bumps, flares and 
plateaus have often been observed in place of smooth power-law decays, thus challenging our 
understanding of the afterglow production. 



7. Conclusions 

We started by testing whether the PLUTO code is appropriate for treating the evolution 
of hyperrelativistic shock waves with Lorentz factors even in excess of 10^, obtainin g satis- 



fact ory evidence of c ohere ncy with self-similar theory developed by iPerna fc Vietril (120021 ) 



and iPalma &: Vietril (120061 ) for shock acceleration in exponential atmosphere. 



To follow we tried to answer the question on how nonlinear effects may let the in- 
stability evolution depart from the linear behavior. We studied several perturbing agent 
configurations, concluding that the shock will tend to restore the original planar shape of 
the discontinuity surface on a time scale given by Tsm ~ 25r/c~^. We intend to remark that 
such a behavior does not appear as a typical saturation phenomenon, due to the lack of a 
competition between a destabilizing factor (which is actually missing) and a restoring agent. 
The reason is easily found: while in the Newtonian counterpart of the problem the desta- 
bilizing factor is given by the tendency of the zero-th order solution to preserve any speed 
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differe nce between adjacent flow columns (even better, it grows indefinitely; see IChevalier 



( ll990l )). here the acceleration is in terms of an homogeneous growth of the Lorentz F factor. 
As a consequence, in the hyperrelativistic regime, essentially due to the existence of the 
speed limit c, even in absence of restoring effects, the maximum gap that a difference of F 
between two distinct cylinders can produce is 

This explains the lack of a substantial destabilizing factor: even if two regions of the shock 
- for example as a consequence of an inhomogeneous upstream - travel with different speeds 
and are at different positions, as soon as the perturbing agent disappears, they will tend 
to reach sooner or later a maximum gap. At that point the only acting process is the 
smoothening influence of the secondary shocks. As a result, the shock will tend inexorably 
to the zero-th order solution (in other terms, the saturation point is that of no perturbation!). 

In the last section we found that self-similarity is reached almost immediately in the 
flow and concluded that this allows us to predict the shock position as a function of the 
only initial and final upstream density. We applied such a result to the case of a turbulent 
ambient density and derived an analytical expression for the dispersion of the shock positions 
at different transverse positions. 

The subject we treated in this paper is expected to have a great relevance with regard 
to models describing GRB radiation, especially those concerning the afterglow emission 
at the external forward shock. The recent Swift observations have show n lots of bump, 



flares and plateaus in place of smooth power-law decays in the light-cu rves (IFox fc Meszaros 
2006). thus ch a llengi n g our understariding o f the afterglow production. iLazzati et al\ ( 120021 ): 



Heyl Sz Pernal (120031 ): iNakar fc OrenI (12004 ) suggested that such a variability in the light- 
curves may be the result of the presence of density bumps in the upstream medium. The 
upcoming launch of GLAST is likely to provide an even more detailed description of such 
events, thus requiring a more accurate modeling of the underlying physics. It goes without 
saying, therefore, that having a good theory for the dynamics of highly relativistic shock 
waves represents a key point of our capability to properly predict the emission expected 
from an afterglow. 

Future developments of this work will involve the evolution of shocks propagating in 
a magnetized ambient medium. An easy and straightforward extension of the numerical 
setup hitherto developed would dispel the uncertainties greatly affecting at the moment this 
theoretical issue and may provide unambiguous evidences on a hot topic of the GRB theory. 
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